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Abstract. We present the results of three-dimensional simulations of super- 
sonic Euler turbulence with grid resolutions up to 1024 3 points. Our numerical 
experiments describe nonmagnetized driven turbulent flows with an isothermal 
equation of state and an rms Mach number of 6. We demonstrate that the iner- 
tial range scaling properties of turbulence in this strongly compressible regime 
deviate substantially from a Kolmogorov-like behavior previously recovered for 
mildly compressible transonic flows. 



1. Introduction 

Understanding the nature of supersonic turbulence is of fundamental importance 
in astrophysics and in aeronautical engineering. In the interstellar medium, 
highly compressible turbulence is believed to control star formation in dense 
molecular clouds. A whole class of more terrestrial applications deals with the 
drag and stability of projectiles traveling through the air at hypersonic speeds. 

Molecular clouds have a highly inhomogeneous structure and the intensity 
of t heir internal motions corresponds to an rms Mach number of the order of 
20. E arso 3 CL981) has demonstrated that within the range of scales from 0.1 pc 



to 100 pc, the gas density and the velocity dispersion tightly correlate with the 
cloud size. 1 Supported by other independent observational facts indicating scale 
invariance, these relationships are often interpreted in terms of supersonic tur- 
bulence with characteristic Reynolds numbers Re ~ 10 8 . Within a wide range 
of densities around 10 3 cm -3 the gas temperature remains close to ~ 10 K since 
the thermal equilibration time at these densities is shorter than a typical hy- 
drodynamic time scale. Thus, an isothermal equation of state can be used as a 
reasonable approximation. While self-gravity, magnetic fields, chemistry, cool- 
ing and heating, as well as radiative transfer should be ultimately accounted for 
in turbulent models of molecular clouds, we focus here specifically on hydrody- 
namic aspects of the problem. 

Numerical simulations of dec aying supersonic hydrodynam ic turbulence with 
the piecewise parabolic method (IColella fc Woodward! fl98l . PPM) in two di- 
mensions were pioneered by iPassot et alJ ((19881 ) 2 and then followed up with 



1 Scc Kapla rTfe Pronikl J1953D for an earlier version of what is now known as Larson's relations. 
2 See a review on compressible turbulence bv lPouauet et all (1991) for references to earlier works. 
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high resolution 2D and 3D simulations bv iPorter et~all (|l992al lbL Il99i [19981. 
ISvtine et alJ ( 20001 ) compared the results of PPM Euler computations with PPM 
Navier-Stokes results and showed that Euler simulations agree well with the high- 
Re limit attained in the Navier-Stokes models. The convergence in a statistical 
sense as well as the direct comparison of structures in configuration space in- 
dicate the ability of PPM to a ccurately sim u late t urbulent flows over a wide 
range of scales. More recently, IPorter et al.l ((20021 1 discussed measures of in- 
termittency in simulated driven transonic flows at M ach numbers of the order 
unity on grids up to 512 3 points. IPorter et al.l (|l999h review the results of these 
numerical studies focusing on the origin and evolution of turbulent structures 
in physical space as well as on scaling laws for two-point structure functions. 
One of the important results of this fundamental work is the demonstration of 
compatibility of a Kolmogorov-type ( Kolmogorovl Il94ll . K41) spectrum with a 
mild gas compressibility at transonic Mach numbers. 

Since most of the computations discussed above assume a perfect gas equa- 
tion of state with the ratio of specific heats 7 = 7/5 or 5/3 and Mach numbers 
generally below 2, the question remains whether this result will still hold for near 
isothermal conditions and hypersonic Mach numbers characteristic of dense parts 
of star forming molecular clouds where the gas compressibility is much higher. 
What kind of coherent structures should one expect to see within the iner- 
tial range of scales in highly supersonic isothermal turbulence? Do low-order 
statistics of turbulence follow the K41 predictions closely in this regime? How 
intermittent is the turbulence? The interpretation of astronomical data from 
new surveys of cold ISM and dust in the Milky Way by Spitzer and Herschel 
satellites requires more detailed knowledge of these basic properties of supersonic 
turbulence. 

In this paper we report first results from our large-scale numerical simula- 
tions of driven supersonic isothermal turbulence at Mach 6 with PPM and grid 
resolutions up to 1024 3 points. We solve the Euler equations in a periodic box 
of linear size L = 1 with initially uniform density distribution p = 1 and the 
sound speed c = 1. We initialize the simulation on a grid of 512 3 points with the 
random velocity field uo that contains only large-scale power within the range 
of wavenumbers k/k m i n £ [1,2], where k m i n = 2ir, and corresponds to the rms 
Mach number M = 6. 

The same velocity field is then used, with an appropriate normalization, as 
a steady random force (acceleration) to keep the total kinetic energy within the 
box on an approximately constant level during the simulation. The random force 

is isotropic in terms of the total specific kinetic energy per dimension, \Uq x 



u o,y) = \ u o,z/i but its solenoidal (V • Uq = 0) and dilatational (V x Uq = 0) 
components are anisotropic since one of the three directions is dominated by the 
large-scale compressional modes, while the other two are mostly solenoidal. The 
distribution of total specific kinetic energy E = lf u 2 dV = E s + E c between the 
solenoidal E s and dilatational E c components is such that xo = Eq/Eq ~ 0.6. 

The driving field is helical, but the mean helicity is very low: (ho) <C (h^), 

where the helicity h is defined as h = u ■ V x u. In compressible flows with 
an isothermal equation of state the mean helicity (h) is conserved, as in the 
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incom pressible case, since the Ertel's potential vorticity is identically zero (Gaffet 
11985ft . 

We first ran the simulation on a grid of 512 3 points from the initial condi- 
tions up to five dynamical times to stir the gas within the box. The dynamical 
time is defined as td = L/(2M). Then we doubled the resolution and evolved the 
simulation for another 5^ on a grid of 1024 3 points. We allow one dynamical 
time for relaxation at high resolution to reach a statistical steady state after re- 
gridding. The time-average statistics are computed using 170 snapshots evenly 
spaced in time over the final segment of 4i^. We use the full set of 170 snapshots 
to derive the density statistics, since the density field displays a very high degree 
of intermittency. This gives us a very large statistical sample, e.g. ~ 2 10 11 mea- 
surements are available to determine the probability density function (PDF) of 
the gas density. The time-average power spectra discussed below are also based 
on the full data set. The velocity structure functions are derived from a sample 
of 20% of the snapshots covering the same period of 4t<2. The corresponding 
two-point PDFs are built on 2 — 4 x 10 9 pairs per snapshot each, depending on 
the pair separation. 



2. Results 

Time evolution of global variables. The time variations of the rms Mach num- 
ber and of the maximum gas density are shown in Fig. ^i, b. The kinetic energy 
oscillates between 18 and 22, roughly following the Mach number evolution. No- 
tice the highly intermittent bursts of activity in the plot of p ma x(t)- The time 

average enstrophy $7 = ~ J \u>\ 2 dV w 10 5 and the Taylor scale A = v/^f ~ 
0.03 = 30 A, where A is the linear grid spacing and u> = V x u is the vorticity. 
The rms helicity grows by a factor of 7.7 in the initial phase of the simulation 
and then remains roughly constant at a level of 1.2 x 10 3 . The conservation of 
the mean helicity is satisfied reasonably well as {h) is contained within ±2% of 
its rms value during the whole simulation. 

Density PDF. In isothermal turbulence the gas density does not correlate with 
the loc al Mach number. As a result, the density PDF follows a lognormal distri- 
bution JVazauez-Semadem|l994 : iPadoan et all 19*9^ : iPassot &; Vazauez-Semadenil 
1998: lNordlund fc Padoanlll999T l. Fig. QJ: shows our results for the time-average 
density PDF and its best-fit lognormal representation. The density distribution 
is very broad due to the very high degree of compressibility in isothermal su- 
personic conditions. The density contrast ~ 10 6 is orders of magn itude higher 
than in the transonic case at 7 = 1.4 studied bv lPorter et alJ ((2002f ). Notice the 
excellent quality of the fit over more than eight decades in probability and a very 
low noise level in the data. If we express the standard deviation a as a function 
of Mach number as a = ln(l + b 2 M 2 ) we get the b e st-fit value of b « 0.26 that 
is smaller than b ~ 0.5 determined in lPadoan et al. I (Il997f ) for supersonic MHD 



turbulence. The powerful intermittent bursts in p m ax{t) correspond to large 
departures from the time-average PDF caused by head-on collisions of strong 
shocks. These events are usually followed by strong rarefactions that are seen 
as large oscillations in the low density wing of the PDF and also in the density 
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Figure 1. Results for the 1024 3 simulation of driven Mach 6 turbulence 
with the ratio of specific heats 7 = 1.001: (a) the time evolution of the rms 
Mach number; (b) the max density as a function of time; (c) the histogram 
for the gas density and the best-fit lognormal; (d) the density power spectrum 
compensated by k; (e) the velocity power spectrum compensated by fc 2 ; (f) 
the power spectrum of dilatational velocity compensated by k 2 ; (g) the 2nd 
order structure functions and the best-fit power laws for log 10 x S [—1.5, —0.5]; 
(h) the 3rd order structure functions and the best-fit power laws for log 10 x £ 
[-1.5,-0.5]. 
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power spectrum. Intermittency is apparently present in supersonic turbulence. 
Its signature in turbulence statistics will be addressed elsewhere. 

Inertial range scaling. The power spectrum of the gas density shows a short 
straight section with a slope of —1.07 ± 0.01 in the range of scales from 250A 
down to 40A followed by flattening due to a power pileup at higher wavenum- 
bers, see Fig. A similar power-law section, although with a slope of — 1.95± 
0.02, and an excess of power in the same near-dissipation range is also clearly 
seen in the velocity power spectrum, Fig. Apparently, the flattening of 
the spectra in the near-dissipation part of the inertial range is due to the 
so-called bottleneck effect related to a three-dimensional non -local mechanism 
of energy transfer between modes of differing length scales (lFalkovichl ll994 N ). 



The same phenome non has been observed both experimentally and in numeri- 
cal simulations (e.s. IPorter et al J 1 1991 iKaneda et ail 1200.1 : IDobler et ai1l2003l : 
iHaueen Sz BrandenburdfeOO^ . At the resolution of 512 3 , the bottleneck effect 
at high wavenumbers and turbulence forcing at low k leave essentially no room 
for the uncontaminated inertial range in the /c-space, even though the spectrum 
does show a nice power law. For the density power spectrum, we would get 
a slope of —0.90 ± 0.02 at 512 3 that is substantially shallower than —1.07 at 
1024 3 . Thus , spec tral index estimates based on low resolution simulations (e.g. 
iKim fc B,vul f2005) may bear a rather large uncertainty. We also note that the 
time-averaging over many snapshots is essential to get the correct slope for the 
density power spectrum since it exhibits variations on a short (compared to 
td) time-scale. The spectrum tends to get shallower upon collisions of strong 
shocks, when the PDF's high density wing rises above its average lognormal 
representation. 

The time-average velocity power spectra for the solenoidal and dilatational 
components show the inertial range power indices of —1.92 ± 0.02 and —2.02 ± 
0.02, respectively, see Fig.^for the fc 2 -compensated dilatational spectrum. Both 
spectra display flattening due to bottleneck with indices of —1.70 and —1.74 in 
the near-dissipative range. The fraction of energy in dilatational modes quickly 
drops from about 50% at k = k m i n down to 30% at k/k m i n « 50 and then 
returns back to a level of 45% at the Nyquist frequency. 

Overall, the inertial range scaling of the velocity power spectrum in our 
highly compressible turbulence simulations tends to be closer to the Burgers 
scaling with a power index of —2 rather than to the K 41 — 5/3-scaling sugg ested 
for the mildly compressible transonic simulations by IPorter et all ^2002). To 
substantiate this result, we also computed the scaling properties of velocity 
structure functions S p (x) = {\u(r +x) — u(r)\ p ) of the orders p = 1, 2, and 3, 
where the averaging is taken for all positions r and all orientations of x within 
the computational domain. Both longitudinal (u \\ x) and transverse (u _L x) 
structure functions can be well approximated by power laws S p (x) oc x^ p , see 
Fig. ^ and h. The low-order structure functions are less susceptible to the 
bottleneck contamination and might be better suited for deriving the scaling 
exponents ( Dobler et al.l I2003T ). 3 The best-fit 2nd order exponents measured 



3 Note, however, that the bottleneck corrections grow with the order p <|Falkovichlll994l) and 
influence the structure functions in a non- local fashion (IDobler et alJl2003D . 
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for the range of scales between 32A and 256A, Q = 0.952 ± 0.004 and C2" = 
0.977 ± 0.008 are substantially larger than the K41 predicted value of 2/3 and 
agree well with our measured velocity power indices. The first order exponents 

(f = 0.533 ± 0.002 and Ci = 0.550 ± 0.004 are significantly larger than the K41 
index of 1/3. We find that the 3rd order scaling exponents Q = 1.26 ± 0.01 and 
£3 = 1.29 ± 0.01 are also noticeably off from unity predicte d by K41 for the 
incompressible limit and also confirmed bv lPorter et al.l ((2002) for the transonic 
regime. Our measurements for the low order exponents roughly agree w i th the 
estimates Ci" ~ 0.5, C2" ~ 0.9, and C3 ~ 1-3 obtained bv iBoldvrev et alJ (|2002T t 
for numerical simulations of isothermal Mach 10 MHD turbulence at a resolution 
of 500 3 points. 

The large-scale driving force we used in these simulations is not perfectly 
isotropic due to the uneven distribution of power between the solenoidal and 
dilatational modes (perhaps, a typical situation for the interstellar conditions). 
We also use a static driving force that could potentially cause some anomalies 
on time-scales of many dynamical times. However, while strong anisotropics 
can significantly af fect the scaling of high-order moments (|Porter et al.l I2OO2I : 
iMininni et al.ll2006l h the departures from Kolmogorov-like scaling we observe in 
the lower order statistics appear to be too strong to be explained solely as a 
result of the specific properties of our driving. The sensitivity of our result to 
turbulence forcing remains to be verified with future high resolution simulations 
involving a variety of driving options. 



3. Conclusions 

Using high-resolution numerical simulations of nonmagnetic highly compress- 
ible turbulence at rms Mach number of 6, we have demonstrated that scaling 
exponents of low-order statistics deviate substantially from Kolmogorov laws for 
incompressible turbulence. A much higher than 1024 3 resolution is required to 
possibly trace a transition from a steeper supersonic inertial range scaling at 
lower A; to a flatter Kolmogorov-like transonic scaling at higher wavenumbers. 

Acknowledgments. This research was partially supported by a NASA 
ATP grant NNG056601G, by NSF grants AST-0507768 and AST-0607675, and 
by a NRAC allocation MCA098020S. We utilized computing resources provided 
by the San Diego Supercomputer Center. 



References 

Boldyrev, S., Nordlund, A., & Padoan, P. 2002, ApJ, 573, 678 
Colella, P. & Woodward, P. R. 1984, J. Comp. Phys., 54, 174 

Dobler, W., Haugen, N. E., Yousef, T. A., & Brandenburg, A. 2003, Phys. Rev. E, 68, 
026304 1 

Falkovich, G. 1994, Phys. Fluids, 6, 1411 
Gaffet, B. 1985, J. Fluid Mech., 156, 141 

Haugen, N. E. & Brandenburg, A. 2004, Phys. Rev. E, 70, 026405 1 
Kaneda, Y., Ishihara, T., Yokokawa, M., Itakura, K., & Uno, A. 2003, Phys. Fluids, 15, 
L21 

Kaplan, S. A. & Pronik, V. I. 1953, Dokl. Akad. Nauk SSSR, 89, 643 



High Resolution Simulations of Supersonic Turbulence 



7 



Kim, J. & Ryu, D. 2005, ApJL, 630, L45 

Kolmogorov, A. N. 1941, Dokl. Akad. Nauk SSSR, 30, 299 

Larson, R. B. 1981, MNRAS, 194, 809 

Mininni, P. D., Alexakis, A., & Pouquet, A. 2006, ArXiv Physics e-prints, 0602148 
Nordlund, A. K. & Padoan, P. 1999, in Interstellar Turbulence, ed. J. Franco & A. Car- 
raminana, 218 

Padoan, P., Nordlund, A., & Jones, B. J. T. 1997, MNRAS, 288, 145 
Passot, T., Pouquet, A., & Woodward, P. 1988, A&A, 197, 228 
Passot, T. & Vazquez-Semadeni, E. 1998, Phys. Rev. E, 58, 4501 
Porter, D., Pouquet, A., & Woodward, P. 2002, Phys. Rev. E, 66, 026301 
Porter, D. H., Pouquet, A., Sytine, I., & Woodward, P. 1999, Physica A, 263, 263 
Porter, D. H., Pouquet, A., & Woodward, P. R. 1992a, Th. Comp. Fluid Dyn., 4, 13 
— . 1992b, Phys. Rev. Lett., 68, 3156 
— . 1994, Phys. Fluids, 6, 2133 

Porter, D. H., Woodward, P. R., & Pouquet, A. 1998, Phys. Fluids, 10, 237 
Pouquet, A., Passot, T., & Leorat, J. 1991, in IAU Symp. 147: Fragmentation of Molec- 
ular Clouds and Star Formation, ed. E. Falgarone, F. Boulanger, & G. Duvert, 
101-118 

Sytine, I. V., Porter, D. H., Woodward, P. R., Hodson, S. W., & Winkler, K.-H. 2000, 

J. Comp. Phys., 158, 225 
Vazquez-Semadeni, E. 1994, ApJ, 423, 681 



